Below is an example of 20 hypothetical 1ha restoration plots (10,000m2, 100m*100m) distributed at random at three reefs in the Townsville region (Davies Reef, Big Broadhurst Reef, Little Broadhurst Reef). The code takes the Allen Coral atlas geomorphology data (see here for more details), filters all reef polygons larger than 1ha in size, randomly selects 20 polygons across the three reefs, and constructs a 100*100m plot around the centroid of each polygon:

# data packages
library(tidyverse)
library(units)
library(foreach)

# spatial packages
library(sf)

# mapping packages
library(tmap)
library(leaflet)

# create function to get make polygons for plots
set_restoration_plot_centres <- function(input = NULL, width = NULL, length = NULL) {

  # Calculate the coordinates of the rectangular polygon
  x <- sf::st_coordinates(input)[1, 1]
  y <- sf::st_coordinates(input)[1, 2]

  # set parameters
  x_min <- x - (width / 2)
  x_max <- x + (width / 2)
  y_min <- y - (length / 2)
  y_max <- y + (length / 2)

  polygon <- sf::st_polygon(list(rbind(c(x_min, y_min), c(x_min, y_max), c(x_max, y_max), c(x_max, y_min), c(x_min, y_min)))) |>
    sf::st_sfc(crs = 20353)

  return(polygon)
}



habitats <- c("Plateau", "Back Reef Slope", "Reef Slope", "Sheltered Reef Slope", "Inner Reef Flat", "Outer Reef Flat", "Reef Crest")

cairns_geomorphic <- st_read("/Users/rof011/coraldynamics/data/Moore-Arlington-20230906220745/Geomorphic-Map/geomorphic.geojson") %>%
   sf::st_transform(crs = sf::st_crs("EPSG:20353")) %>%
   dplyr::group_by(class) %>%
   dplyr::mutate(habitat_id = paste0(
   gsub(" ", "_", class), "_",
   sprintf(paste0("%0", ceiling(log10(max(1:length(class)))), "d"), 1:length(class)))) %>%
   sf::st_make_valid() %>% 
#  mutate(habitat_id=as.factor(habitat_id)) %>%
  dplyr::group_by(habitat_id, class) %>% 
  summarise(geometry = st_union(geometry)) %>% 
  mutate(area = round(st_area(geometry))) %>%
  filter(class %in% habitats) %>%
  drop_na(class) 
## Reading layer `Moore Arlington' from data source `/Users/rof011/coraldynamics/data/Moore-Arlington-20230906220745/Geomorphic-Map/geomorphic.geojson' using driver `GeoJSON'
## Simple feature collection with 3768 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 145.9182 ymin: -16.89574 xmax: 146.2582 ymax: -16.61326
## Geodetic CRS:  WGS 84
cairns_plots_centroids <- cairns_geomorphic %>%
  filter(area > set_units(10000,"m^2")) %>%
  filter(class %in% c("Reef Slope", "Sheltered Reef Slope", "Back Reef Slope")) %>% 
  drop_na(class) %>%
  st_centroid() %>% 
  st_cast("POINT")


cairns_plot_outputs <- foreach(i=1:nrow(cairns_plots_centroids), .combine="c") %do% {
  tmp <- set_restoration_plot_centres(cairns_plots_centroids$geometry[i], width=100, length=100)
  tmp
}

Map restoration plots

Select layers using the layercontrol on the left, use [ ] for full-screen viewing.

set.seed(1)

cairns_plot_outputs <- cairns_plot_outputs |> st_as_sf() |> mutate(id=paste0("Plot_ID_",seq(1,n(),1))) 

cairns_plot_outputs2 <- cairns_plot_outputs[sample(nrow(cairns_plot_outputs), 20), ] # subset to 20

# visualise with thematic maps
map <- tm_view() +
tm_tiles("Esri.WorldImagery", group = "<b> [Seascape]</b> satellite map", alpha = 0.5) +
  tm_shape(cairns_geomorphic |> mutate(class=as.factor(class)), id="area", name = "<b> [Seascape]</b> habitats") +
  tm_borders(col = "black", lwd = 0.2) +
  tm_fill("class", name = "Benthic habitats", id="area", palette = c("Plateau" = "lightskyblue1", "Back Reef Slope" = "darkcyan",
                                                                    "Reef Slope" = "darkseagreen4", "Sheltered Reef Slope" = "darkslategrey",
                                                                    "Inner Reef Flat" = "darkgoldenrod4", "Outer Reef Flat" = "darkgoldenrod2",
                                                                    "Reef Crest" = "coral3"), alpha = 0.6) +
  tm_shape(cairns_plot_outputs2, id="plot_id", name = "<b> [Restoration]</b> plots") +
    tm_fill(col="darkred", alpha=0.6) +
    tm_borders(col="black")


  map |> tmap::tmap_leaflet() |>
    leaflet::addProviderTiles('Esri.WorldImagery', group = "<b> [Seascape]</b> satellite map", options=leaflet::providerTileOptions(maxNativeZoom=18,maxZoom=100)) |>
    leaflet::addProviderTiles('Esri.WorldTopoMap',  group = "<b> [Seascape]</b> base map", options=leaflet::providerTileOptions(maxNativeZoom=19,maxZoom=100)) |>
    leaflet::addLayersControl(position="topleft", overlayGroups=c(
      "<b> [Seascape]</b> satellite map", "<b> [Seascape]</b> habitats", "<b> [Restoration]</b> plots"),
                              options=leaflet::layersControlOptions(collapsed = FALSE)) |>
    leaflet.extras::addFullscreenControl(position = "topleft", pseudoFullscreen = FALSE)

Restoration area

What does 20ha cover across the three reefs? Using the spatial data we can calculate restoration area as a proportion of reef habitat (combined reefs: combined reefs: Moore, Thetford, Arlington) assuming 20ha is located across i) Total coral area, ii) Total reef slope area, iii) Total reef flat area:

Habitat Area Proportion
Total coral area 85598086 0.2%
Total reef slope area 28713698 0.7%
Total reef flat area 56884388 0.4%